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Comment on "Mechanical analog of temperature for the description of force 
distribution in static granular packings" 

Philip T. MetzgeiEl 

The KSC Applied Physics Laboratory, The John F. Kennedy Space Center, NASA 
YA-C3-E, Kennedy Space Center, Florida 32899 
(Dated: February 2, 2008) 

It has been proposed by Ngan [Phys. Rev. E 68, 011301 (2003)] that the granular contact 
force distribution may be analytically derived by minimizing the analog of a thermodynamic free 
energy, in this case consisting of the total potential energy stored in the compressed contacts minus 
a particular form of entropy weighted by a parameter. The parameter is identified as a mechanical 
temperature. I argue that the particular form of entropy cannot be correct and as a result the 
proposed method produces increasingly errant results for increasing grain rigidity. This trend is 
evidenced in Ngan's published results and in other numerical simulations and experiments. 

PACS numbers: 45.05.+X, 45.70.-n, 81.05.Rm, 05.70.-a 



Ngan recently proposed a functional minimization 
approach to derive the granular contact force probabil- 
ity distribution. The method is different from the en- 
tropy maximization approaches that have been proposed 
by others because it explicitly accounts for the poten- 
tial energy stored in the compressed granular contacts. 
The previous methods had assumed Shannon's entropy 
for either the distribution of contact force Cartesian 
components fx (where a; is a principal stress axis) , 

/•OO 

S = -k dU PAU) lnP.(/.), (1) 
Jo 

or the distribution P of contact force magnitudes / Q , 

/"OO 

S=-k d/P(/) lnP(/). (2) 
Jo 

In either case, the entropy was maximized subject to the 
conservation of the number of contacts and the conserva- 
tion of stress(es) in the packing. This, along with other 
important variations of the methods, predicted a distri- 
bution function. It is the latter of these two entropies 
that Ngan adopts. In a paper now in preparation jjj I 
show why both of these entropies are incorrect, because 
the density of states in the relevant phase space must 
be profoundly non-uniform, and it is this non-uniformity 
which is the source of the unique shape of the contact 
force distribution function in the region of weak forces. 
In this Comment I will discuss the problem with the en- 
tropy as it is relevant to Ngan's hypothesis and show how 
the problem is revealed in the published results. 

The functional proposed by Ngan is F = U — OS where 



U 



d/ P{f) W{f) 



(3) 
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is the internal energy of the packing and where W{f) 
is the work function for compression of a contact under 
normal force /. Thus, U is an analog of the Helmholtz 
free energy from thermodynamics and the parameter 9 is 
the proposed mechanical analog of temperature. This 
hypothesis produces interesting results because, for a 
Hertzian contact force law, the predicted distribution 
function closely fits the numerical simulation data both 
for 2D and, in the case of strong deformations, for 3D (at 
least to the statistical precision of the numerical data). 
My own work has focused on a theoretical analysis of the 
case of perfect rigidity which is the limit in which Ngan's 
method is incorrect. I am not able to comment on the 
region of large deformations, which I suspect is where 
Ngan's work makes its key contribution. 



I. GAUSSIAN VERSUS EXPONENTIAL 

In Ngan's hypothesis the form of the tail depends upon 
W{f), predicting for Hertzian contacts a Gaussian tail 
^ eicp{—Pf'^) in 2D or a nearly Gaussian compressed ex- 
ponential tail ~ exp(— /3/^/'^) in 3D. This is in substantial 
agreement with statements in a paper by O'Hern, Langer, 
Liu, and Nagel (OLLN) 5j, which Ngan cites. Before ad- 
dressing the form of the entropy in Ngan's hypothesis, it 
is necessary to question whether Gaussian tails have re- 
ally been observed in static granular force distributions. 
I claim that they have not, except for cases with very 
small numbers of grains or large deformations 0, 13- 
Five arguments make this case. 

First, the tail normally does not appear exponential 
(a straight line on a semilogarithmic plot) until several 
multiples of the average value of force, (/). This onset 
is apparently dependent upon the dimensionality of the 
packing. For example. Figs. 2 and 5 of Ref. 9 show the 
onset at / « 3 (/) for 2D and / « 2 (/) for 3D. There 
may be other factors which hasten or delay the onset as 
well. A number of the semilog plots in the published 
literature lack statistical precision for a sufficient range 
beyond the onset and thus may appear Gaussian. This is 



2 



because the eye tends to extrapolate the curvature it sees 
just prior to the loss of precision. However, I have seen 
no case that is truly inconsistent with an exponential tail 
beyond the reasonable range of onsets except when there 
are large deformations or a very small number of 

grains ^Gj. On the other hand, there are numerous exam- 
ples which clearly show an exponential tail inconsistent 
with a Gaussian or other curved form. These examples 
include 2D frictional simulations using contact dynamics 
(CD) and molecular dynamics 3D frictional sim- 

ulations using CD 13, a Hertzian contact law and 
a Hookcan contact law [Til IT^ ; 3D frictionless simula- 
tions using two versions of the Lennard-Jones potential 
[T^j l: 3D frictional experiments 0|; and 3D frictionless 
experiments with emulsions • 

Second, the paper by OLLN claimed that 3D friction- 
less static packings with a harmonic potential produce a 
Gaussian tail j^. The thermal argument predicting this 
tail was not correctly applied to static packings. The 
density of states of a thermal ensemble and the density 
of states of a static packing ensemble are organized by 
two mutually exclusive principles: one by the conserva- 
tion of energy and momentum with respect to time, the 
other by the conservation of forces with respect to several 
spatial dimensions. In the transition as the temperature 
T — > 0, a packing does not become static until it locates 
one of the relatively rare locations in phase space that 
satisfy the highly organizing static equilibrium require- 
ments. There is no basis to assume that the density of 
states, or the potential energy distribution representing 
it, will have the same form after undergoing this self- 
organizing transition. 

Third, the empirical data from OLLN's paper do 
not contradict the existence of an exponential tail, either. 
Their two Figs. (3) (a) and (3)(b) represent two types of 
ensemble averages. In the first, an exponential tail onset- 
ting at / « 3 (/) fits arguably better than OLLN's pro- 
posed Gaussian over the same region In the second, 
the tail is clearly exponential with an onset at / ~ 2 (/). 
For the harmonic potential in which / = Kx, renormaliz- 
ing the force scale of a packing is equivalent to renormal- 
izing its spatial scale. Hence, the first ensemble can be 
interpreted as all possible packing geometries compressed 
to achieve the same average force at the expense of dif- 
ferent packing fractions. The second ensemble is the set 
of all possible packing geometries compressed to achieve 
the same packing fraction but at the expense of different 
average forces. Arguably it is the latter ensemble aver- 
age, not the former, which represents the self-averaging of 
a much larger packing. This is because locally averaged 
force fluctuations do occur in large packings (because the 
spatial distribution of forces is dominated by force chains 
and is therefore very heterogeneous), whereas the mean 
offset in OLLN's Fig. (2)(b) implies that locally averaged 
packing fraction fluctuations become relatively small in 
an increasingly large packing. Hence, this implies that an 
exponential onset in very large, frictionless, 2D packings 
ought to occur closer to / « 2 (/) , which is in agreement 



with the other data cited above. However, regardless of 
which ensemble average is the "correct" one, both are at 
least consistent with an exponential tail. I believe this 
conclusion is in better agreement with the statements 
found in the follow-on full-length paper by O'Hern, Sil- 
bert, Liu, and Nagel p^ . 

Fourth, a simple argument can prove that, for fric- 
tionless packings in the limit of small deformations, the 
force distribution cannot be a functional of W{f). Be- 
cause a frictionless packing is isostatic flTj , all the forces 
can be derived deterministically from the geometry of the 
contact network and the imposed boundary loads, alone. 
Therefore, the only role that W{f) can play is by help- 
ing to decide what the geometry of the contact network 
shall be. If W{f) is relatively soft so that deformations 
are large, then the geometry of the contact network will 
be perturbed signiflcantly, and then W{f) may indeed 
affect the resulting forces. However, in the limit when 
the deformations become vanishingly small, then the ge- 
ometry of the contact network becomes independent of 
the form of W{f) and the contact forces also become 
independent of the form of W{f). Therefore, the force 
distributions resulting from every W{f) must approach 
a universal form in this limit. 

Fifth, Bouchaud's argument with regard to the 
q model persuasively explains the universal form of the 
tail. He shows that the sufficient condition for an expo- 
nential tail is merely that some grains be allowed to tip 
all of their loads into one contact, which allows arbitrarily 
large forces to accumulate along particular force chains. 
This argument applies to all cohesionless granular me- 
dia regardless of W{f), except for cases where the grains 
cannot freely tip. We may draw an important inference: 
when deformations become large, the formation of addi- 
tional contacts increases the stability of the grains, hence 
allowing the forces to be more evenly distributed through 
space _8j. This may erode the exponential tail so that it 
transits to a more rapidly decreasing form |19| . 



II. THE ROLE OF THE ENTROPY 

These arguments do not imply that Ngan's free-energy 
hypothesis is necessarily incorrect. W{f) is a measure 
of deformation, which is the relevant parameter when 
the tail becomes nonexponential. The hypothesis does 
produce excellent results in the case of large deforma- 
tions. The tail of the predicted distribution also becomes 
steeper and more curved as the deformations increase, 
thus matching the observed trend X^]. However, it can 
be seen in Ngan's 3D results, Fig. (7), that the analyti- 
cal predictions do not fit the simulation data for the 3D 
case with the least hydrostatic pressure applied to the 
packing, i.e., the case with the least grain deformation. 
In fact, there is an unmistakable trend that with smaller 
pressures the predicted tail is increasingly distant from 
the simulation data which become increasingly consistent 
with an exponential tail (2l|. It is possible with just a 
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small change in Ngan's hypothesis to make the predicted 
tail follow the trend toward an exponential. This could 
be done by modifying the definition of the mechanical 
temperature 6 so that k remains constant and 9 oc 
Then, the influence of W{f) would vanish where the iso- 
static argument says that it should, in the limit where 
deformations are small. 

However, as shown in Ngan's Figs. (1) and (2), 
straightening the tail would raise a secondary but more 
fundamental problem because it would produce the incor- 
rect features in the region of weak forces. Then, because 
these features affect the average value of the distribution, 
the straightened tail would have the wrong exponential 
decay constant, (3 = (/) instead of /? « 1.6(/). Like 
the exponential tail, these weak-force features are so uni- 
versal that they, too, are most likely the result of some 
fundamental, organizational behavior in granular media. 
They are the small peak near the average value of force 
(or a plateau, or at least an abrupt change in slope in 
cases of severely anisotropic stress |23) and the nonzero 
probability density at zero force. 

This secondary problem can be related to the Shan- 
non's entropy maximization hypothesis of Kruyt and 
Rothenburg (KR) 3]. If Ngan's U ^ Q with a nonva- 
nishing fc as I have suggested, then the functional F in 
Eq. ^ reduces to the entropy term alone. Minimizing 
this is KR's method but with frictional grains. If we 
took KR's prediction in the limit as the Coulomb coeffi- 
cient of friction vanishes, — > 0, then it, too, becomes a 
purely exponential distribution, contrary to published re- 
sults in the region of weak forces for frictionless packings 

The problem lies with Shannon's entropy because of 
its inherent implication that all possible sets of contact 
forces are equally probable, meaning that the density 
of states in the phase space with coordinates {/i | i ~ 
1, . . . , is uniformly populated over the entire accessi- 
ble region. The accessible region is implied by the use of 
the Lagrange multipliers to be all states where the aver- 
age force per contact is correct, or in the case of Ref.|3 all 
states where the volumetrically averaged stresses match 
the externally applied stress tensor. To be more accurate, 
the Shannon's entropy implies that any non-uniformity 
in the accessible region will be unbiased with respect to 
the distribution of coordinates, so that the weaker forces 
are not left out by the nonuniformities any more often 
than are the stronger forces. In granular packings this as- 
sumption is fundamentally incorrect because it neglects 
the most important organizing feature of granular pack- 
ings: the requirement that the grains be stable. There 
is no analogous requirement in classical thermal systems, 
and it turns out that this requirement biases the density 
of states against weaker forces. 

An explanation of the bias begins by noting that any 
two contact forces on the same grain are strongly cor- 
related, increasingly so as the contacts are further away 
from each other toward opposite sides of the grain 
For simplicity this Comment must illustrate how the cor- 




FIG. 1: Special case to illustrate why grains with below- 
average forces correspond to fewer stable locations in phase 
space than do grains with average forces. 

relation affects the density of states using only a spe- 
cial case, allowing the reader to draw the connections 
to the general case. For an isotropic packing of fric- 
tionless 2D grains, consider one grain which has two of 
its contacts exactly opposite one another, ipi — Q and 
'!/'2 = Ti" as shown in Fig. ^ The force /i is clearly re- 
lated to /a, but not generally equal to it because of the 
contact forces /2 and fi located in 7r/3 < -01 < 27r/3 
and 47r/3 < -04 < 57r/3, respectively. Note that steric 
exclusion keeps the contacts arranged fairly predictably 
around a grain, which prevents the statistically averaged 
general case from deviating too far from this special case. 
Static equilibrium requires 

/2 sin 02 = -/4sini/'4, 

h ^ h~ G, (4) 

where 

G = /2 cos 02 + fi cos 04 (5) 

is the difference between /i and /s. Then G = J~^f4, 
7-1 / / ^^"^^4 

J = cos V'4 — COS ^02 ;— • (6) 

sm 02 

Neglecting that this is a special case (for illustrative sim- 
plicity), we assume that the distribution of is repre- 
sentative of all the forces in the packing. These have a 
distribution Pf{f), which we want to derive. By chang- 
ing variables the distribution of G in terms of Pf is 

Pg{G)= / d02 / dV'4 |J| P/(G J), (7) 

where J is identified as the Jacobian. All we know 
about Pf is that it is zero for negative arguments (tensile 
forces), but positive valued and normalized for positive 
arguments. Assuming it is a continuous function its nor- 
malization implies that it has a bounded tail. Over the 
range of integration, J has odd symmetry in the sense 
that J (-02, "04) = — J (tt — -02, 37r — 04). Hence, Pg{G) 
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must be an even function which is positive valued over 
— oo < G < oo, having a bounded tail in both extremes. 
We cannot solve the distribution of /i by directly chang- 
ing variables from (G, /s) to (/i, f^) because the pairs are 
not statistically independent and we do not know their 
joint probability distributions. However, we know that 
the only stable configurations of this grain are those in 
which G < /a because of Eq. Q . We can therefore cal- 
culate the proportions of the volume of phase space that 
include only the stable configurations of this grain with 
a particular value of f^. Since there are two equations 
of stability there are only two degrees of freedom in the 
force dimensions. Integrating across only one of these to 
find the size of stable space as a function of the other, 



v{h) « yydV d/4 e{h)Pf{h) 

oc /" dG e(/3-G)PG(G) 

'J —OO 

oc / dG Pg (G) , 

J — OO 



(8) 



where 9 is the Heaviside(unit step) function. Hence, V 
is a monotonically increasing function of /a which has a 
finite value at /a = 0. This illustrates that the volume of 
the stable regions in phase space is biased against weaker 
forces, but is not vanishing for zero force. Neglecting 
that this is a special case, the Pf corresponding to the 
maximum volume of stable phase space is, therefore. 



Pf{f)=e 



-Pf 



f 



which is a recursion equation in Pf through Eq. JT)). 
Thus, the relative slope of the exponential factor and the 
integral will determine the behavior of Pf in the region of 
weak forces, and this provides the mechanism to explain 
variations in that region as a function of the anisotropy 
of the packing [2^ . 



dG Pg(G), 



(9) 



Because of the bias against weak forces in stable phase 
space, the Shannon's entropy adopted by Ngan and oth- 
ers is not correct and this explains my claim why Ngan's 
predicted form for the case of vanishing deformations 
cannot fit both the tail and the region of weak forces 
for any choice of 9 and k. If we could define the en- 
tropy so that it accommodates the nonuniform density 
of states, then perhaps a particular choice of 9 and k can 
fit both the tail and the region of weak forces even in the 
low-deformation limit. However, in the case of very de- 
formable grains where Ngan's hypothesis works best, the 
principal role of the entropy term seems to be that it in- 
jects a logarithm into the equations. Whether or not we 
will ever deduce a correct entropy maximization method 
for the case of perfect rigidity, it is safe to say that the 
actual form of the entropy with the relevant Lagrange 
multipliers will still contain a logarithm, so the applica- 
bility of Ngan's work in the large deformation limit is 
still an open and interesting question. 

I am grateful to A.H.W. Ngan for the helpful interac- 
tion which has significantly improved this Comment. 
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